Evaporation of quark drops during the cosmological quark— hadron transition 
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We have carried out a study of the hydrodynamics of disconnected quark regions during the final 
stages of the cosmological quark-hadron transition. A set of relativistic Lagrangian equations is 
presented for following the evaporation of a single quark drop and results from the numerical solution 
of this are discussed. A self-similar solution is shown to exist and the formation of baryon number 
density inhomogeneities at the end of the drop contraction is discussed. 
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I. INTRODUCTION 

According to the standard picture of cosmology, a 
phase transition during which a plasma of free quarks 
and gluons was converted into hadrons should have oc- 
curred at about 10 -5 s after the Big Bang. At present 
it is not completely clear whether this transition would 
have been first order or of a higher order [1], although 
recent lattice calculations favour it being first order for 
the realistic case of two degenerate light u and d quarks 
and a heavier s quark (of up to 400 MeV) [2]. The in- 
teresting consequences that a first order transition could 
have produced for the early history of the universe [3,4] 
provide the motivation for studying this case. 

The work described here is within the scenario of 
a first-order transition starting with the nucleation of 
hadronic bubbles in a slightly supercooled quark-gluon 
plasma. The bubbles are thermodynamically favoured 
and grow with the phase interface moving as a deflagra- 
tion front and with the quark-gluon plasma being pro- 
gressively transformed into a plasma of hadrons. The 
dynamical time scale of the growth is likely to be much 
shorter than the universe expansion time scale, and so 
the hadronic bubbles grow until they meet and eventu- 
ally coalesce. This next stage of the transition is proba- 
bly characterized by disordered motions and dissipative 
processes which transform the kinetic energy associated 
with the ordered motion of the quark regions into inter- 
nal energy of the two phases. The temperature jump 
between the two phases becomes extremely small with 
both phases being almost at the critical temperature T c 
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and the rate of transformation of the quark-gluon plasma 
being controlled by the overall expansion of the universe. 
When more than half of the space is occupied by hadronic 
matter, regions containing quark matter start to become 
disconnected and, while evaporating, they tend to be- 
come spherical under the effects of surface tension. 

We concentrate here on the hydrodynamics of a single 
evaporating quark drop, making use of the mathemati- 
cal and numerical experience developed for the previous 
study of growing hadronic bubbles [5-8]. A particular 
feature of our study is a demonstration of the existence 
of a self-similar solution which is likely to be relevant for 
some of the final stages in the shrinking away of quark 
drops and which provides the basis for our discussion of 
the final stages of the transition. These stages are par- 
ticularly important since it is then that relics could be 
produced which might then survive until later epochs. 
The possible astrophysical consequences of a first order 
transition have been extensively investigated in the liter- 
ature both with regard to the formation of baryon num- 
ber inhomogeneities [9-11] and their possible effects on 
primordial nucleosynthesis and also for the production 
of primordial magnetic fields [12,4]. We will comment 
on these expectations and discuss the implications of the 
results obtained here for discussions of the formation of 
baryon number inhomogeneities. 

The following is a brief outline of the contents of the 
paper. The set of relativistic hydrodynamical equations 
and junction conditions at the deflagration front are re- 
viewed in Section II. The existence and derivation of the 
self-similar solution for a shrinking drop is discussed in 
Section III, while Section IV is dedicated to the presen- 
tation of the method for the computation. Here we com- 
ment on the specification of the initial conditions and 
discuss numerical tests which have been carried out. Re- 
sults from the computations are presented in Section V, 
where we also discuss their consequences for the forma- 
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tion of baryon number density inhomogeneities. Finally, 
conclusions and perspectives for future developments of 
the work are presented in Section VI. In this paper, we 
use units for which c = k = k = 1 and denote standard 
partial derivatives with a comma. 



II. RELATIVISTIC HYDRODYNAMICAL 
EQUATIONS 



radiative transfer. This was described in [5,6] to which 
the reader is referred for further details. 

Consider a Lagrangian formulation in which the space- 
time metric has the spherically symmetric line element 



ds 2 



-a 2 dt 2 + b 2 dp 2 + R 2 (d6 2 + sin 2 6d<p 2 ) , (1) 



with p being a comoving radial coordinate having its ori- 
gin at the centre of symmetry. The system of hydrody- 
namical equations then can be written as 



We here consider a single, isolated, spherical quark 
drop, whose surface is inward-moving as a consequence of 
the evaporation. The drop has initial dimensions which 
are much smaller that the typical distance between dis- 
connected quark regions (often taken to be of the order 
of 10 14 — 10 16 fm), but larger than the shortest mean free 
path of the particles which are not strongly-interacting. 
These are electromagnetically and weakly interacting 
particles (mainly photons, muons, electrons, neutrinos 
and their respective antiparticles), and although they 
are not directly involved in the phase transition, they 
can provide a long-range energy and momentum trans- 
fer at a scale comparable with their mean free path A. 
For neutrinos and antineutrinos A ~ 10 13 fm, while for 
electromagnetically interacting particles A ~ 10 4 fm. A 
quark drop with initial dimensions of the order of 10 5 
fm is essentially transparent to the neutrinos (which will 
be neglected here), while the strongly interacting mat- 
ter is in thermal and mechanical equilibrium with the 
electromagnetically interacting particles. Because all of 
these are essentially massless particles, we can consider 
them as components of a generalized "radiation fluid" 
and their dynamical and energetic contribution is taken 
into account by making a suitable modification of the 
number of degrees of freedom in the equations of state 
for the quark and hadron phases. 

The thermal and mechanical equilibrium ceases to hold 
when the drop has reached dimensions comparable with 
the mean free path of the radiation fluid which then be- 
comes progressively more decoupled from the strongly- 
interacting matter. Although a complete treatment of 
the evaporation of a quark drop, involving the solution of 
the relativistic radiative transfer problem together with 
the pure hydrodynamical one, is in principle possible, 
many of the main features of the evaporation process al- 
ready emerge in a computation in which the radiation 
fluid is always considered as completely coupled to the 
strongly interacting ones. Performing a calculation in 
which long-range energy and momentum transfer is taken 
into account is certainly a more difficult task and making 
an error analysis without first having a clear underlying 
hydrodynamical solution would be much more compli- 
cated. For these reasons the set of hydrodynamical equa- 
tions used here is not the one presented in [7,8], where 
the relativistic radiative transfer problem was solved for 
a growing bubble, but rather the equivalent one without 



AitR'-p^ 
w 



in 

G ( — + 4irpR 



Rt 
P R 2 


= au, 
- a U)li 

R,n ' 
= w P,t, 

_ e )|U - wp tfl 


aw 


pw 




= AitR 2 eR )ll , 


r 


= 4tt P R 2 R^ : 


b 


= (4ttR 2 p)~ 1 



.1/2 



(2) 
(3) 

(4) 
(5) 

(6) 

(7) 
(8) 
(9) 



Here u is the radial component of the fluid four- velocity in 
the associated Eulerian frame, R is the Schwarzschild cir- 
cumference coordinate, T is the general relativistic ana- 
logue of the Lorentz factor, p is the relative compression 
factor, w is the specific enthalpy (w = (e + p)/p, with e 
being the energy density and p the pressure). For small 
net baryon number, it is reasonable to take e and p as 
depending only on temperature, and then the equations 
of state for the two phases of the strongly interacting 
matter can be written as 
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where the hadron medium is considered as consisting of 
massless, point-like pions, while the quark phase is de- 
scribed by the Bag Model equation of state, with B being 
the "Bag" constant. We here take g q = 37, gu = 3 and 
g e = 9, where g e accounts for the degrees of freedom of 
the electromagnetically interacting particles. 

As in the hadronic bubble growth computation, we 
treat the transition region as a discontinuity surface with 
the fluid and metric variables on either side being linked 
by suitable junction conditions. Making use of the Gauss- 
Codazzi equations [13], the junction conditions for energy 
and momentum can be expressed respectively as 
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(12) phase interface is to be treated as a discontinuity, then 
it seems that the only satisfactory way of accomplish- 
ing this is by using a system of characteristic equations 
[v_ into which the hydrodynamical equations (2)-(9) can be 

ab transformed. For this we have used the equations and 

methodology presented in [6], to which the reader is re- 
; (13) ferred for further details. 



where p, s = dp: s /dt is the interface velocity (ji s < 0), 
/ = (a 2 — b 2 /! 2 ) 1 ! 2 and a is the surface tension. The 
labels ± indicate quantities immediately ahead of the 
transition front (quark phase) and immediately behind it 
(hadron phase) and [A] ± = A + - A~ , {A}* 1 = A + + A~ . 
Other supplementary junction conditions follow from 
continuity across the interface of the metric quantities 
R, dR/dt, and ds 
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and from the time evolution of the jump in the mass 
function m 
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For the circumstances under consideration here, it is sat- 
isfactory to calculate the jump of the mass function at 
the interface as [m] ± = 4irR 2 tj. Because the drop sur- 
face moves as a deflagration front, one further equation 
is needed expressing the rate at which the quark matter 
is transformed into hadrons at the phase interface. A 
suitable expression is obtained by setting the hydrody- 
namical flux F H into the hadron region equal to the net 
thermal flux F„ into it 
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Here a is an accommodation coefficient (0 < a < 1) 
which can be associated with the "transparency" of the 
phase interface to the thermal flux and is, at least in prin- 
ciple, calculable from theory. Hereafter we will assume 
a = 1. A discussion of the behaviour of solutions for 
lower values of a can be found in [6]. (Note that in de- 
riving (18) we have taken the unit space-like four-vector 
normal to the interface as the one pointing towards the 
centre of the drop.) 

A main requirement for the numerical solution of the 
system of hydrodynamical equations and junction con- 
ditions at the interface is that there should be a correct 
causal treatment of the deflagration front [6,7]. If the 



III. SELF-SIMILAR SOLUTIONS 

In general a self-similar motion can be expected to oc- 
cur when there are no intrinsic length or time scales in 
the system. When this is the case, the hydrodynamical 
equations can be re-written in terms of a single dimen- 
sionless variable obtained as a suitable combination of 
the parameters of the problem. 

Self-similar growth of bubbles in first order cosmologi- 
cal transitions has been considered in the literature both 
in the case of detonation fronts [14-16] and deflagration 
fronts [17,18]. For bubble expansion, self-similar motion 
appears only when the bubble radius is small enough so 
that there is no interaction between bubbles, but large 
enough so that surface tension effects are negligible. In 
the case of drop contraction we expect that self-similar 
motion should set in when the drop radius is reasonably 
smaller than the mean distance between quark regions, so 
that one can neglect the interaction between drops, but 
large enough so that surface tension effects are negligible. 

In this section, it is convenient to work with the special 
relativistic form of the hydrodynamical equations writ- 
ten in terms of Eulerian coordinates (here denoted by 
r and t). Gravity does not play any important role on 
the scale of the drops being considered here and so it 
is quite sufficient to consider only the special relativis- 
tic equations for discussing these self-similar solutions. 
For spherical symmetry, the standard conservation and 
continuity equations are 

A [(e+^ 2 ) 7 2 ] +A [(e+^T 2 ^] 
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where j = (1 — i; 2 ) -1 / 2 is the relativistic Lorentz factor 
and v is the three velocity of the fluid. In the case of a 
one-parameter equation of state, such as the one which 
we are using, it is sufficient to solve just equations (19) 
and (20) and equation (21) is then needed only if one 
wishes to calculate p. For a relativistic system, the only 
possible choices for a similarity variable are £ = ± r/t 
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[19]. We take £ = r/t and write e = e(£), p = /?(£) and 
v = «(£). In terms of £ the system (19)— (21) becomes: 



where c s = (dp/ de) 1 ! 2 = (1/3) 1 / 2 is the sound speed. 
The similarity variable £ = r/t can be viewed either as 
the position of a point in the solution at a given time, or 
as the velocity at which a given point in the profile of the 
solution is moving. This velocity is to be clearly distin- 
guished from the fluid velocity v at the point described 
by £. A self-similar flow which is linearly expanding with 
time is naturally described in terms of a positive simi- 
larity variable £. Conversely, a self-similar contracting 
flow, such as the one occurring for an evaporating quark 
drop, is described in terms of negative values of £. In this 
case, the time can be thought of as progressing through 
negative values and tending to zero as the radius of the 
contracting drop tends to zero. 

When solving equations (22)-(24) numerically, it is 
useful to notice their property of invariance under the si- 
multaneous transformations £ — ► — £ and v —> —v which 
has the consequence that it is only necessary to solve 
them in one of the half-planes £ £ [0,± 1] in order to 
know the solution in the whole interval £ £ [— 1, 1]. In 
Figs. 1 and 2, we have plotted the results of numerical 
integrations for the functions v(£) and e(£). The dashed 
lines in Fig. 1 represent points for which the derivative 
dv/d£ in (22) becomes infinite, i.e. for which 

(cte - iy + 2*(i - c > + ( C ? - e) = o, (25) 

which has the roots 



Expression (26), which is also the special relativistic for- 
mula for velocity composition, expresses the fact that 
the fluid velocity, measured relative to an observer mov- 
ing at £, equals the speed of sound. A physical solution 
for v cannot be extended across these points, since it 
would then be double valued, and a discontinuity must 
be introduced. The upper right quadrant of Fig. 1 shows 
the solutions of the similarity equations for an expanding 
system (such as a spherically expanding hadron bubble 
[17,6]). The upper left quadrant shows the equivalent so- 
lutions for a contracting system (such as an evaporating 
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FIG. 1. Solutions of the similarity equations for the veloc- 
ity v. The relevant quadrant for a contracting drop is the 
upper left-hand one. The point £ = —c s represents the sonic 
radius, while the centre of the drop is at £ = at any given 
time. Parts of solution curves similar to the dot-dashed one 
are used in construction of the similarity solutions shown in 
Fig. 3, where they appear in a mirror image. 
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FIG. 2. Solutions of the similarity equations for the the en- 
ergy density e. The solutions have been suitably normalized 
to the value of the energy density at the sonic point and cor- 
respond to solutions in the upper quadrants of Fig. 1. The 
relevant curves for a contracting drop are in the negative part 
of the £ plane. The point £ = —c s represents the sonic ra- 
dius, while the centre of the drop is at £ = at any given 
time. Parts of solution curves similar to the dot-dashed one 
are used in construction of the similarity solutions shown in 
Fig. 4, where they appear in a mirror image. 
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spherical quark drop), while the lower quadrants pro- 
vide the corresponding solutions for negative values of 
the fluid velocity (note that these do not necessarily rep- 
resent physically realistic configurations). 

For the contracting quark drops which we are consid- 
ering here, the point £ = represents the centre of the 
drop at any given time, £ = — c s is the "sonic radius" and 
£ = — 1 is the edge of the past light cone. The only self- 
similar flow solution in Figs. 1 and 2 which can be taken 
to extend for all £ is the trivial solution v = 0, e = const. 
Any other solution must be the result of a "patching" of 
different solution curves with the joins being either via a 
weak discontinuity (where the function is continuous but 
has a derivative of some order which is not continuous), 
or via a full discontinuity (where the function itself is not 
continuous). The first applies in the case of the edge of a 
rarefaction or compression wave, while the second occurs 
in the case of shocks or combustion fronts. 

The situation for a contracting drop is radically dif- 
ferent from that for an expanding bubble. For the con- 
tracting drop, the deflagration front is moving inwards 
and spherical symmetry imposes that the only possible 
state ahead of it is one with zero velocity and constant 
density since the only self-similar solution satisfying the 
condition v = at £ = is the trivial one. At the drop 
surface, a discontinuity corresponding to the deflagration 
front must be introduced in order to join onto the relevant 
self-similar solution curve for the flow behind it. Once Th 
has been specified at the interface (i.e. T^), the junction 
conditions give the velocity of the front and the state of 
the fluid adjacent to it, which correspond to points in the 
(v, £) and (e, £) planes of Figs. 1 and 2. The solution for 
the medium behind the front then follows the continuous 
curves passing through the relevant points. 

For a weak deflagration (the physically interesting 
case), the self-similar solution behind the front is charac- 
terized by decreasing velocity, which goes to zero at the 
sound radius (£ = — c s ), and increasing energy density. 
At the sound radius, the solution then joins (via a weak 
discontinuity) onto another one with zero velocity and 
constant energy density. Note that a weak discontinu- 
ity, (which propagates at the local sound speed) can only 
be located at the point £ = — c s and the sonic radius is 
then the edge of the perturbed flow region at any given 
time. These weak deflagration solutions are illustrated 
with the solid line curves in Figs. 3 and 4. Different 
curves are drawn for different values of the temperature 
in the hadron phase adjacent to the phase interface, with 
= /T c ) ranging between 0.95 for the leftmost solid 
curve and 0.61 for the rightmost solid curve. The latter 
represents a Chapman-Jouguet deflagration, for which 
the front moves at the sound speed relative to the fluid 
behind [20]). This marks the transition from weak de- 
flagrations (which are physically realistic) to strong de- 
flagrations (which are physically forbidden and are in- 
dicated with dashed lines in Figs. 3 and 4). Note that 
the rightmost solid curve and the leftmost dashed curve 
correspond to adjacent but distinct values of . 
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FIG. 3. Self-similar curves for the velocity for different val- 
ues of the temperature in the hadron phase adjacent to the in- 
terface. The solid curves represent weak deflagration solutions 
with the rightmost solid curve being a Chapman-Jouguet de- 
flagration. The dashed curves represent strong deflagration 
solutions beyond Chapman-Jouguet limit and are physically 
forbidden. 
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FIG. 4. Self-similar curves for the energy density for dif- 
ferent values of the temperature in the hadron phase adja- 
cent to the interface. The solid curves represent weak de- 
flagration solutions with the rightmost solid curve being a 
Chapman-Jouguet deflagration. The dashed curves represent 
strong deflagration solutions beyond Chapman-Jouguet limit 
and are physically forbidden. 
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The solution for the medium behind a weak deflagra- 
tion front can be seen as a compression wave since the en- 
ergy density is increasing to the value in the unperturbed 
fluid. This occurs because a fluid element crossing the 
phase interface experiences a strong decompression from 
its original state and is subsequently compressed to the 
background value. Conversely, for the strong deflagra- 
tion solutions, the medium behind the front is acceler- 
ated away, reaching a velocity approaching c at infinity, 
while the energy density decreases through a rarefaction 
wave, going to zero at infinity. Note that as decreases 
from 1, the temperature of the quark phase first decreases 
but then reaches a minimum (at = 0.70) and starts 
to increase again. This is consistent with the fact that 
the function T q = T q (Th) is not monotonic (see Fig. 3 of 

When analyzing self-similar deflagration solutions for 
contracting drops, a number of novel features are encoun- 
tered. Three main differences with respect to the equiv- 
alent solutions for expanding bubbles are: 

i) In the case of an expanding bubble, the deflagration 
is always preceded by a compression wave fronted, 
in principle, by a shock (which has negligible ampli- 
tude in the case of a hadronic bubble nucleated with 
small supercooling [6,8]). For a contracting drop, 
the deflagration can only be preceded by a solution 
with zero velocity and constant energy density. 

ii) The solution for the medium behind the deflagra- 
tion front in an expanding bubble could differ from 
a constant state with zero velocity only if the bub- 
ble were expanding supersonically, in which case 
it would have to be followed by a rarefaction wave 
[18]. The medium behind the deflagration in a con- 
tracting drop is never at rest, but rather is outflow- 
ing with the velocity tending either to zero at the 
sonic point (in the case of weak deflagrations) or to 
the light velocity at infinity (in the case of strong 
deflagrations). 

Hi) Depending on the degree of supercooling, the de- 
flagration front in an expanding bubble could move 
with velocities up to the sound speed (and possi- 
bly beyond) relative to the centre of symmetry. Ir- 
respective of whether the deflagration is weak or 
strong, the surface of a contracting drop always 
moves very subsonically relative to the centre of 
symmetry (i.e. at the surface, |£| <C c s ). 

The existence and properties of the self-similar solu- 
tions for a contracting drop will be used when specifying 
initial conditions for the numerical solution of the full 
hydrodynamical equations (as described in the following 
Sections). The initial conditions used are similar to the 
leftmost solid curves in Figs. 3 and 4. 



IV. NUMERICAL STRATEGY 

A. Organization of the grid 

For following the evaporation of quark-gluon drops, 
we have used a numerical code constructed on the ba- 
sis of the strategies adopted previously for calculating 
the growth of hadronic bubbles [5-7]. As in the previous 
computations, the present code makes use of a compos- 
ite numerical technique in which a standard Lagrangian 
finite-difference method is used to solve the hydrody- 
namical equations in the bulk of both phases, while a 
system of characteristic equations and a set of junction 
conditions are solved in the grid zones adjacent to the 
phase interface. We have used a spherically symmetric 
Lagrangian grid having comoving coordinate /j, and with 
its origin at the centre of the drop. The grid has variable 
spacing with the width of each zone being twice that of 
the zone inside it (i.e. A/j,j + i/2 = A^j.^), apart from 
the two central zones which have equal width. This ar- 
rangement allows bubble expansion or drop contraction 
to be followed through changes of many orders of magni- 
tude. 

A major difference with respect to the previous cal- 
culations is that we are here following an inward-moving 
phase interface rather than an outward-moving one. This 
has necessitated a different organization for the solution 
of the system of characteristic equations and the use 
of modified algorithms for the calculation of quantities 
in the grid zones immediately adjacent to the interface. 
Moreover, in order to keep constant the number of grid 
zones within the drop (as required for maintaining accu- 
racy), it has been necessary to implement a "re-gridding" 
routine which creates a new zone inside the quark phase 
every time the interface crosses a zone boundary during 
its inward motion. 

Regridding can be quite a delicate matter and, if the 
implementation is not a good one, there is a danger of 
introducing instabilities into the solution, particularly if 
function fitting routines are used. With our choice of grid 
structure, it is possible to avoid the use of any fitting al- 
gorithm. Our procedure is that every time the interface 
crosses a zone boundary (in the inward direction) the 
central zone is divided into two equal parts and all of the 
other zones are re-labelled (i.e. j — ► j + 1 ). This strat- 
egy maintains intact the structure of the grid (the first 
two central zones still have the same width, while the 
others are increasingly spaced) and limits the recalcula- 
tion of new quantities only to the central zones. To keep 
constant the total number of zones (and thus avoid in- 
creases in the computational time), the outermost zone 
is removed every time the central zone is divided into 
two. Experience has shown that the loss of information 
involved in this does not create any problem since the 
outer edge of the grid is always very distant from the 
interface. 

This re-gridding scheme has several advantages: it is 
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extremely simple to implement, it is suitable for use dur- 
ing changes of many orders of magnitude in the drop 
radius, and it keeps constant the number of zones within 
the drop, hence maintaining resolution and numerical ac- 
curacy. 



B. Initial Conditions and Numerical Tests 

The specification of initial conditions in a time depen- 
dent calculation can be considered as a "problem within 
the problem" . When there is no obvious physical solu- 
tion, it can be useful to extract information from the 
code itself by observing its response to trial initial data. 
This has been the case with the present calculation, in 
which the initial conditions (which have turned out to 
correspond to self-similar solutions), were originally "sug- 
gested" by the code itself. 

Before initial data based on a self-similar solution was 
implemented, other choices of initial data were observed 
to relax, after some readjustments, to a solution having 
a flat density profile and zero velocity within the quark 
region, and an outward velocity profile falling off roughly 
as R~ 2 in the hadron region. The subsequent evolution of 
this solution exhibited self-similar features and this then 
stimulated the investigation for a similarity solution. It 
is important to stress that it could happen that a self- 
similar solution which is mathematically possible might 
nevertheless not be realized in practice. In this respect, 
the weak deflagration similarity solution has turned out 
to be strikingly "robust". As a test, initial conditions 
with an inward-pointing velocity field inside the drop 
and zero velocity outside, were specified on the zero-time 
hypersurface. After some excursions, the solution con- 
verged towards the corresponding similarity one. Identi- 
cal results have been obtained also in the case when the 
solution was started with very irregular and noisy initial 
conditions. 

During the intermediate stages of the transition the 
universe has been re-heated by the latent heat released 
by the quark confinement, and the temperatures in both 
phases have become close to the critical temperature for 
the transition. The disordered motions resulting from 
the bubble percolation and coalescence tend to produce 
homogeneity in the temperature profile in both phases 
and the dynamics of the transition is then driven by the 
universe expansion, whose cooling effect is compensated 
by the latent heat released. During the last stages of 
the transition, the reduced quark volume fraction is no 
longer able to provide the amount of energy necessary to 
keep the increased hadron volume fraction at T c . As a 
consequence Th is free to decrease and a small non-zero 
temperature jump drives the following hydrodynamical 
evolution of the disconnected quark regions. In view of 
this, we assume that all quantities have reached a high de- 
gree of homogeneity in each phase and then the situation 
is regulated by the temperature jump between the two 



phases. All of the variables, except for the fluid velocity 
u, are therefore taken to have "step-like" profiles. Note 
that this is consistent with the self-similar solution for 
small supercooling and low interface velocity (see Fig. 4). 

Working within this scenario, we first specify the drop 
radius i? s and the temperature of the hadron phase Th 
and the corresponding temperature in the quark phase 
T q is then calculated from the special relativistic form of 
the junction conditions (see eq. (15) in [21]). The next 
quantity to calculate is the value of the interface velocity 
jj, s which is obtained from (18) as 



Ps 



-Ai:R s p + a + 



where 



X 



aTr 2 (g h + g e ) 



(i + x 2 ) 1/2 -i 



(T*-Tt) 



60 



(e+p + 47r 2 ffe T 4 /90) H 



(27) 



(28) 



and, using the freedom in the choice of reference values 
for p and a, we have set p + = a + = 1. Rewriting the 
junction condition (12) as 



P- 



(e+p)- 

P+ , , N X , 
( e +P) + 



(29) 



where x = a_/a + and combining this with (16), it is 
possible to obtain the following equation 



h +Ps 



(e+p)+ fb+p s 



(e+p). 



(30) 



whose solution provides the required value for a_ and 
consequently for p_ . The last quantity to calculate is 
m_, and the expression for this follows directly from (15) 



Ps 



(h-b-), 



(31) 



where, following the similarity solution, we have taken 
u+ = and also u = for all points inside the quark 
drop. From the similarity solutions for weak deflagrations 
shown in Fig. 3, it can be seen that the velocity profile 
just behind the phase interface can be well approximated 
by a solenoidal flow, for which 



uR 



const = m_R 



(32) 



This approximation is an excellent one in the case of 
small supercooling. Expression (32) is applied out to the 
point where the value for u given by this becomes less 
than that for the background universe, which is taken to 
follow the spatially flat Robertson-Walker solution. The 
expression for u then follows from (8), with T being taken 
equal to 1, and is 
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f2Gr 



V R 



1/2 



(33) 



Equations (27)-(33) provide all of the data required on 
the zero-time hypersurface. 

For the present consideration of drop evaporation at 
the cosmological quark-hadron transition, it is expected 
that the degree of supercooling would be rather small 
and this case is our primary interest. However, in or- 
der to gain a good overall perspective of the situation, it 
is interesting to consider also what would happen if the 
supercooling were greater and so we have investigated 
values of Th ranging between 0.999 and the Chapman- 
Jouguet point (Th = 0.61). We found it convenient and 
satisfactory to continue to specify the initial data as out- 
lined above also in these cases. 



V. RESULTS 

When studying the hydrodynamics of an evaporating 
quark drop it can be useful to compare this process with 
the evaporation of a water drop. The main difference be- 
tween the two situations is that while the first is exother- 
mic, the second is endothermic. This means that the wa- 
ter drop will continue to evaporate as long as there is an 
efficient transport of energy from the inner regions to the 
surface where it will be lost during the evaporation. In a 
similar way, it is important that during the quark drop 
evaporation there should be an efficient mechanism in 
the hadron phase for transporting the latent heat away 
from the phase interface. The geometry of the process 
and the complexity of the phenomena at the deflagration 
front, prevent any simple but accurate analysis being pos- 
sible and only a numerical simulation can give a deeper 
insight into the details of the evolution. The present 
calculation is the first complete relativistic computation 
of an evaporating cosmological quark drop, although a 
more simplified analysis was carried out by Kajantie and 
Kurki-Suonio [22]. Numerical integrations of the hydro- 
dynamical equations discussed in the previous sections 
over the whole permitted range of values for the initial 
Th , seem to indicate a scenario for the final stages of 
quark drop evaporation which is rather different from 
the one proposed in [22]. 

We will first present the results obtained for the fol- 
lowing combination of values for the free parameters in 
our model: we consider an initial quark drop of radius 
i? s = 10 5 fm, surrounded by a hadron plasma at temper- 
ature Th = 0.99, to which is associated a phase interface 
with surface tension parameter <7o = <r/T^ = 1. (We 
take T c = 150 MeV.) Results with other initial values 
for R s and Th show broadly similar hydrodynamical be- 
haviour and the choice for the value of Th is related to the 
fact that the computational times increase greatly for Th 
closer to 1. We first discuss the situation for <j = 1 
as this corresponds to the value used in our previous 



T h /T c =0.99 

(7o=1.00 



1 2 3 4 5 6 

logio R (fm) 

FIG. 5. Time evolution of the profile of the energy density 
e. The hadron phase is on the right of the vertical disconti- 
nuity. Note that the similarity solution is preserved until to 
values of the drop radius of the order of 10 2 fm. 
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FIG. 6. Time evolution of u, the radial component of the 
fluid four-velocity in the Eulerian frame. 
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computations [6,8] and also allows one to see clearly the 
main features of the hydrodynamical scenario. Results 
for lower values of <7o will be discussed later. 

During the initial stages of the drop evaporation, the 
evolution of the hydrodynamical variables in both phases 
is in very good agreement with the similarity solution. As 
the evaporation proceeds and the drop reaches dimen- 
sions of the order of R s k, 10 2 fm, the hydrodynamics 
starts to deviate significantly from the similarity solution 
(the evolution is no longer scale-free). The surface ten- 
sion starts to cause the compression inside the shrinking 
drop to increase, raising the temperature jump between 
the two phases which has been constant up to this stage. 
The value oiT^ is not significantly changed by this, how- 
ever. This deviation from the similarity solution is at the 
origin of a run-away mechanism which amplifies the tem- 
perature difference between the two sides of the interface 
and increases the quark evaporation rate. Although the 
area of the drop surface is reduced, the velocity of the 
evaporating matter has become larger and this preserves 
a significant outward flux away from the surface. As a 
consequence, the contraction is accelerated and the drop 
experiences an increasingly rapid evaporation which ends 
with its complete disappearance (see Figs. 5-7) 

In this scenario then, the final stages of the evapora- 
tion are not regulated by the expansion of the universe, 
as suggested by [22], but rather by the run-away mecha- 
nism which allows the temperature jump across the phase 
interface to increase. One way of understanding this pro- 
cess is by taking into account the larger outward flow pro- 
duced by a higher evaporation rate from the drop surface. 
This can, in turn, increase the efficiency in the removal 
of the latent heat released, preventing the temperature 
in the hadron phase rising above T c , while allowing the 
temperature in the quark phase to grow above T c . The 
hydrodynamical behaviour described above has been ob- 
served for all of the allowed values of the hadron temper- 
ature and although the use of lower initial values for Th 
(to which correspond larger temperature jumps between 
the two phases) produces a more dramatic and rapid evo- 
lution for the shrinking quark drop, all of the results ob- 
tained share the above features. A physical limit to the 
values of Th is given by Th = 0.61 which corresponds to 
a deflagration front moving at the sound speed with re- 
spect to the medium behind (i.e. a Chapman- Jouguet 
deflagration). 

The deviation of the hydrodynamical solution away 
from the self-similar one can be seen as a consequence 
of the existence of a surface tension a. When the drop 
dimensions become comparable with the characteristic 
length / rs c/ie-q ~\~Pq)j the drop evolution can no longer 
be considered as scale-free and surface effects start to 
become dynamically relevant. When <7o is taken to be 
exactly zero, the similarity solution continues to hold all 
the way down to the complete disappearance of the drop. 
A numerical computation of this situation has been per- 
formed and the self-similar solution was preserved with 
a precision close to the sixth decimal place. If <7q = 0.02, 
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FIG. 7. Time evolution of the compression factor p. Note 
that during the contraction of the drop from 10 fm to 1 fm, a 
depression appears in the hadron phase. This depression will 
later be compensated by a rarefaction wave. 
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TABLE I. The relative change (AQ/Q ) of the hydrody- 
namical variables for different values of the surface tension (To . 
The relative variation is evaluated between the final and the 
initial values of the relevant variable which is shown in the 
first column. Here m s is the radial component of the inter- 
face four-velocity as measured by an Eulerian observer. All 
variations have been calculated for an initial Th = 0.99. 
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as suggested by recent lattice-gauge simulations [23], the 
breaking of the similarity solution starts to become rel- 
evant later, for drop dimensions R s ps 10 fm, with the 
subsequent evolution following the lines described for the 
case Co = 1. A main difference between these two situ- 
ations is given by the magnitude of the relative change 
in the hydrodynamical variables as the contraction pro- 
ceeds. Bearing in mind that the sooner the self-similar 
solution is broken, the more the variables can grow away 
from their initial values, we report in Table 1 the total 
relative change of the most relevant quantities as calcu- 
lated between the initial value (i.e. the value at the time 
the computation is started) and the final value (obtained 
when R s ps 1 fm) , for different values of <7o . A graphical 
representation of the behaviour of the radial component 
of the interface four- velocity in the Eulerian frame u s , 
and of T q h is also given in Figs. 8 and 9. We note that 
the approximation of treating the phase interface as an 
exact discontinuity certainly breaks down before the drop 
radius reaches 1 fm and so our results for the smallest 
radii should only be treated as indicative. 

One way of looking at the results in Table 1 is that 
of considering the self-similar solution as a perfect bal- 
ance between competing effects. As long as the simi- 
larity holds, there is a dynamical equilibrium by means 
of which the evaporation flux reduces the dimensions of 
the quark drop without increasing the compression or the 
energy density inside it. The closest analogue of this pro- 
cess could be that of a "leaky" balloon which, although 
shrinking, maintains the same shape and internal com- 
pression by means of ejecting internal material. When 
the self-similar solution is broken, the balance is lost and 
the run-away mechanism sets in. Although these very 
final stages in the life of the quark drop represent only 
a small fraction of its whole evolution, they can be ex- 
tremely important since a considerable change in the fluid 
variables might occur then. 

A comment should be made on the results obtained 
for low initial values of the hadron temperature (i.e. 
Th < 0.80). Although these cases probably have no 
cosmological relevance, as they require a strongly su- 
percooled quark plasma, they provide important infor- 
mation about the hydrodynamics of weak deflagration 
fronts near the stability threshold given by the Chapman- 
Jouguet point. In the present scenario this occurs for an 
initial value Th = 0.61 and has turned out to be the 
limit for the numerical solution. While a detailed anal- 
ysis would go beyond the scope of this paper, it is in- 
teresting to note, when looking at the sequence of the 
different solutions in the range of low Th , that the hy- 
drodynamical evolution becomes less stable and regular 
as values of Th closer to the Chapman-Jouguet point are 
used. It is particularly interesting to note the behaviour 
of u s , which is larger for lower initial values of Th . In this 
case, the presence of oscillating modes is very clear and 
these can perhaps be related to the instabilities which 
are known to appear from the linear perturbation anal- 
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FIG. 8. Behaviour of m s for different values of the surface 
tension parameter ao. Note that the deviation away from 
the self-similar solution becomes extremely small for smaller 
values of ao- 
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FIG. 9. Behaviour of the temperature in the two phases 
Th,q for different values of the surface tension parameter ao- 
Note that the deviation away from the self-similar solution 
becomes extremely small for smaller values of ao- 
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ysis of two-dimensional deflagration fronts [24,25]. At 
the Chapman- Jouguet point, the solution is extremely 
unstable and is almost immediately dominated by insta- 
bility modes with very large amplitude and very short 
dynamical time scale. As stated before, it was not possi- 
ble to calculate any numerical solution beyond this limit. 
Figs. 10 and 11 show the behaviour of u s for Th = 0.80 
and Th = 0.61. It is important to point out that, despite 
the dramatic oscillatory behaviour seen for u s at low Th , 
the solution for the other quantities remains compara- 
tively smooth and regular up until the Chapman-Jouguet 
point, beyond which stability is lost completely. This sort 
of behaviour seems to be characteristic of deflagrations 
and is related to phenomena observed in laboratory ex- 
periments [26]. The complete loss of stability past the 
Chapman-Jouguet point (i.e. for strong deflagrations), 
may be related to the fact that in this case the front 
ceases to be in mutual causal connection with the fluid 
behind (relative to which the front is supersonic). 

Much of the astrophysical interest in a first order 
quark-hadron phase transition concerns baryon num- 
ber segregation within the shrinking quark drop and the 
consequent generation of inhomogeneities in the baryon 
number density. The occurrence of baryon number seg- 
regation is related to the fact that in chemical equilib- 
rium, baryon number in the high temperature phase is 
carried by almost massless quarks, while in the low tem- 
perature phase it is carried by nucleons, whose rest mass 
m T c and whose number density is limited by a fac- 
tor exp(— m/T c ). The degree of inhomogeneity might 
then be further increased during evaporation of the quark 
drops. Several studies have been carried out (see [27] for 
a review) in order to calculate the amplitude of the in- 
homogeneities which could be produced, but the lack of 
a detailed knowledge of the micro-physical processes oc- 
curring at the interface (where suppression mechanisms 
could, in principle, reduce the flow of baryon number 
from the high temperature phase), and the difficulty in 
estimating the typical scale at which the hadronic bub- 
bles meet, make these calculations extremely tentative. 

With the present calculation, in which we neglect the 
effects of long-range energy and momentum transfer, we 
have aimed to establish whether, by means of purely hy- 
drodynamical mechanisms, the contraction could signifi- 
cantly aid baryon concentration. With regard to this, a 
key result is the demonstration of the existence of self- 
similar solutions which tend to keep constant the values 
of the compression factor in the two phases. If one makes 
the (very idealised) assumption that baryon number is 
carried along exactly with the hydrodynamical flow and 
limits attention to material which was within the discon- 
nected quark regions at the time when chemical equilib- 
rium was broken, then the baryon number density will be 
found to follow the same behaviour as the compression 
factor. As can be seen from Fig. 7, no strong enhance- 
ment of the compression factor seems to appear. For a 
standard quark drop of <7q = 1 and initial Th = 0.99, 
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FIG. 10. Time evolution of the radial component of the 
interface four-velocity in the Eulerian frame for large su- 
percooling. This should be compared with Fig. 8, where 
T h /T c = 0.99. 
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FIG. 11. Time evolution of the radial component of the 
interface four-velocity in the Eulerian frame for large su- 
percooling. This should be compared with Fig. 8, where 
T h /T c = 0.99. 
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FIG. 12. Space-time diagram of the final stages of the drop 
evaporation. Here i? SOJV marks the sonic radius, while the 
thick solid line is the world-line of the phase interface. The 
oblique straight lines emanating from the point of disappear- 
ance of the drop, represent the rarefaction wave (r.w.) and 
the lines marked with arrows are fluid flowlines. 

the relative increase of the compression at the end of the 
contraction is Ap/ p ~ 0.49 (see Table 1). Slightly larger 
compressions have been found for lower values of Th (e.g. 
Ap/ p = 0.55 for Th = 0.80), but such low values of the 
temperature seem unlikely to be relevant in the cosmo- 
logical situation. 

Our conclusion is that hydrodynamical compression on 
its own can do little to enhance baryon number concen- 
tration. If there is to be a concentration of the magnitude 
discussed in the literature [11,4], then this must result 
mainly from the role of the long-range energy transfer 
or from suppression of the flux of baryon number across 
the interface. However, it is worth noticing that the final 
disappearance of the quark drop will be followed by the 
emission of a rarefaction wave as shown schematically in 
the space-time diagram in Fig. 12. This will act in the 
direction of reducing any overdensity which could have 
formed and its effect needs to be properly taken into ac- 
count. A more extended treatment is now in progress 
involving the solution of the relativistic radiative trans- 
fer problem together with the solution of the deflagration 
hydrodynamics for a contracting drop. (The equivalent 
calculation was carried out in [7,8] for an expanding bub- 
ble.) 

VI. CONCLUSION 

We have presented the relativistic hydrodynamics of a 
contracting quark drop during the cosmological quark- 
hadron phase transition in the absence of long-range en- 
ergy and momentum transfer. We have shown the ex- 



istence of a self-similar solution which governs the hy- 
drodynamics of a contracting quark region in the ab- 
sence of intrinsic length scales. Numerical solution of 
the hydrodynamical equations for a representative range 
of values within the allowed parameter space, has led 
to a new and consistent scenario for the final stages of 
the phase transition in the absence of long-range energy 
and momentum transfer. An isolated quark drop with 
an initial radius large enough so that surface effects can 
be neglected, has been numerically followed during the 
contraction. Initially the evaporation is governed by the 
similarity solution which has the peculiarity of preserving 
the values of the energy density and of the compression 
factors in both phases. When the bubble radius becomes 
sufficiently small (i.e. R s ps c/(e q +Pq)), surface effects 
start to become relevant and the self-similar solution is 
then broken. The drop then experiences an increasingly 
rapid evaporation, eventually ending with the emission 
of a rarefaction wave when the drop finally disappears. 
This hydrodynamical behaviour has been observed for 
values of the temperature in the hadron phase ranging 
between Th = 0.999 and Th = 0.61, which is the sta- 
bility threshold for a weak deflagration front (i.e. the 
Chapman- Jouguet point). 

The simulations performed seem to indicate that in 
the absence of long-range energy and momentum trans- 
fer and with baryon number being entirely carried along 
with the hydrodynamical flux, no relevant baryon num- 
ber concentration appears above that indicated from con- 
siderations of chemical equilibrium in the middle part of 
the transition. This study provides a hydrodynamical 
background for a more complete treatment in which long- 
range energy and momentum transfer via electromagnet- 
ically interacting particles will be taken into account. 
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